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We  connect  explicitly  the  classical  0(2)  model  in  1  +  1  dimensions,  a  model  sharing  important  features 
with  U(l)  lattice  gauge  theory,  to  physical  models  potentially  implementable  on  optical  lattices  and  evolving 
at  physical  time.  Using  the  tensor  renormalization-group  formulation,  we  take  the  time  continuum  limit  and 
check  that  finite-dimensional  projections  used  in  recent  proposals  for  quantum  simulators  provide  controllable 
approximations  of  the  original  model.  We  propose  two-species  Bose-Hubbard  models  corresponding  to  these 
finite-dimensional  projections  at  strong  coupling  and  discuss  their  possible  implementations  on  optical  lattices 
using  a  S7Rb  and  41 K  Bose-Bose  mixture. 
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I.  INTRODUCTION 

Recently,  there  has  been  a  lot  of  interest  in  the  possibility 
of  building  quantum  simulators  for  lattice  gauge  theory  (LGT) 
using  optical  lattices  [1-5].  The  goal  is  to  engineer  many-body 
systems  with  cold  atoms  that  can  be  built  experimentally  and 
that  approximately  evolve  according  to  some  given  quantum 
LGT  Hamiltonian.  Achieving  this  goal  would  allow  us  to  go 
beyond  what  can  be  done  with  classical  computing,  namely 
overcoming  the  sign  problem  of  quantum  chromodynamics 
(QCD)  with  a  chemical  potential,  establishing  its  phase 
diagram,  and  studying  its  real-time  evolution.  Introducing  a 
chemical  potential  in  QCD  is  necessary  to  describe  physical 
situations  in  which  a  nonzero  quark  density  is  needed,  such  as 
the  early  universe  or  heavy-ion  collisions.  Building  a  quantum 
simulator  for  QCD  requires  that  we  first  systematically 
establish  the  viability  of  the  approach  by  building  up  on  simple 
models  sharing  some  of  the  basic  features  of  lattice  QCD. 

In  the  context  of  condensed  matter,  a  proof  of  principle  that 
quantum  simulating  is  possible  has  been  given  in  the  case  of 
the  Bose-Hubbard  model.  For  this  simple  model,  a  remarkable 
level  of  quantitative  agreement  [6]  has  been  reached  between 
state-of-the-art  quantum  Monte  Carlo  calculations  and  their 
experimental  optical  lattice  implementations.  It  would  be  very 
desirable  to  provide  a  similar  proof  of  principle  in  the  context 
of  LGT. 

In  this  article,  we  propose  an  optical  lattice  setup  and 
accurate  numerical  methods  to  relate  it  to  a  simple  model 
that  shares  some  important  features  (discrete  imaginary  time, 
relativistic  space-time  symmetry,  compact  gauge  variables, 
and  a  complex  action)  with  interesting  LGT  models,  namely 
the  classical  0(2)  in  1  +  1  dimensions  with  a  chemical 
potential.  This  model  is  described  in  Sec.  II.  The  goal  of  the 
article  is  to  discuss  the  optical  lattice  implementation  of  one  of 
the  building  blocks  of  the  Hamiltonian  formulation  of  gauge 
theory,  namely  the  “quantum  rotors”  that  are  described  in 
more  detail  below,  rather  than  discussing  more  specific  aspects 
such  as  the  implementation  of  Gauss’s  law  for  LGT  models 
involving  these  building  blocks. 


The  connection  between  the  classical  0(2)  in  1  +  1  dimen¬ 
sions  and  physical  systems  on  optical  lattices  requires  three 
steps.  First,  we  introduce  computational  methods  based  on 
the  tensor  renormalization-group  (TRG)  method  [7-10]  to  take 
the  time  continuum  limit  (step  1,  Sec.  II)  and  to  calculate 
the  effects  of  finite-dimensional  truncations  necessary  for  a 
physical  implementation  (step  2,  Sec.  III).  We  then  construct 
a  two  species  Bose-Hubbard  model  which  at  second  order  in 
degenerate  perturbation  theory  can  be  matched  with  the  finite¬ 
dimensional  truncations,  and  we  propose  an  experimental 
implementation  using  a  87  Rb  and  41 K  Bose-Bose  mixture 
(step  3,  Sec.  IV).  The  0(2)  model  is  very  well  understood 
using  classical  computing  [7-12],  and  our  goal  is  not  to  learn 
more  about  this  model  from  quantum  simulations  but  rather  to 
demonstrate  that  a  quantitative  correspondence  is  possible. 

One  should  be  aware  of  the  fact  that  in  contrast  to  the 
quantum  Monte  Carlo  treatment  of  condensed-matter  models 
where  space  and  time  are  completely  independent  entities, 
the  state-of-the-art  calculations  in  LGT  are  performed  using 
the  Lagrangian  formalism  at  discrete  imaginary  time  where 
space  and  time  are  completely  interchangeable.  In  LGT,  the 
continuum  limit  is  usually  taken  in  a  way  that  preserves 
this  relativistic  symmetry  between  space  and  time.  The 
Hamiltonian  representation  provides  the  functional  forms  used 
to  fit  correlation  functions,  and  a  slightly  better  resolution  in  the 
time  direction  is  sometimes  used,  however  the  time  continuum 
limit  is  not  taken  independently.  Explicit  Hilbert  space  repre¬ 
sentations  of  the  physical  states  and  of  their  matrix  elements 
are  mostly  absent  from  today’s  lattice  QCD  calculations.  In  our 
construction,  the  first  step  will  be  to  take  the  time  continuum 
limit  using  the  Lagrangian  formulation.  Note  that  Lorentz 
symmetry  can  emerge  near  criticality  in  the  Hamiltonian 
formulation  [13]  and  that  the  classical  0(2)  model  is  often 
used  as  an  effective  theory  for  the  Bose-Hubbard  model  [14]. 

It  is  important  to  understand  the  similarity  between  the 
infinite-dimensional  Hilbert  spaces  of  the  0(2)  model  and 
U(l)  LGT  in  the  Hamiltonian  formulation.  In  the  mid-1970s, 
LGTs  were  developed  in  the  Hamiltonian  formalism  [15-18] 
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using  local  gauge  variables  that  live  on  bonds  connecting 
neighboring  sites.  For  continuous  and  compact  symmetry 
groups,  these  gauge  links  are  operators  that  live  on  an 
infinite  Hilbert  space  and  in  the  appropriate  basis  look  like 
classical  group  elements.  In  U(l)  LGT  [16],  gauge  links  are 
phases  e  ,  which,  when  considered  as  operators,  live  in  an 
infinite-dimensional  Hilbert  space  spanned  by  the  eigenstates 
| n)  of  the  “angular  momentum”  operator  L  —  —id/86  with 
all  positive  and  negative  integer  eigenvalues  n.  The  same 
“quantum  rotors”  appear  in  the  Hamiltonian  formulation  of 
the  0(2)  model  [17,18], 

For  realistic  implementations  with  cold  atoms,  it  is  neces¬ 
sary  to  consider  Hamiltonians  where  gauge  links  are  quantum 
operators  that  live  in  a  finite  rather  than  infinite  Hilbert 
space  [19,20].  In  the  U(l)  example,  this  would  mean  the 
eigenvalues  of  L  only  take  a  finite  range  of  values.  For  this  to 
occur  naturally,  one  restricts  the  Hilbert  space  to  be  in  a  spin- 
s  representation,  i.e.,  n  =  —s,  — (s  —  1), ...  0,  ...  (s  —  l),s. 
Finite-dimensional  projections  and  quantum  link  variables 
have  played  an  important  role  in  recent  proposals  to  simulate 
dynamical  gauge  fields  [1-3,21,22], 

The  common  features  of  the  0(2)  model  considered  here 
and  the  U(l)  gauge  model  can  be  understood  by  comparing 
the  TRG  formulations  of  the  two  models  [8],  In  both  cases, 
the  Fourier  expansion  of  exp  [/3  cos(<9)]  is  used,  which  leads 
to  the  labeling  of  states  by  (positive  and  negative)  integers. 
However,  the  quantum  numbers  are  associated  with  plaquettes 
in  the  gauge  case  rather  than  links  in  the  spin  case.  The  physics 
of  the  models  is  also  quite  different.  For  instance,  in  2  +  1 
dimensions,  the  0(2)  spin  model  has  a  second-order  phase 
transition  while  the  U(l)  gauge  model  has  none. 

n.  THE  MODEL  AND  ITS  TIME  CONTINUUM  LIMIT 

The  simplest  model  involving  the  quantum  rotors  described 
above  is  the  0(2)  model  in  1  +  1  dimensions.  Its  partition 
function  reads 


with  action 

S  —  ^  '  cos(6(VJ  j  ]  lO 

(*■0 

-  ps  ^2  cos(6>(jc+1,o  -  0(x, o)-  (2) 

(x.t) 

The  meaning  of  the  chemical  potential  //  [23]  appears  clearly 
in  the  limit  where  f3s  is  zero  and  we  have  decoupled  quantum 
rotors  with  a  discrete  spectrum  labeled  by  nx  at  each  site  x  [see 
Eq.  (5)].  Using  these  labels,  the  chemical  potential  generates  a 
contribution  —finx  to  the  energy  at  each  site.  The  sites  of  the 
rectangular  Ns  x  Nr  lattice  are  labeled  as  ( x,t )  and  we  assume 
periodic  boundary  conditions  in  space  and  time. 

When  )$>  ps,  we  obtain  the  time  continuum 
limit  [17,1 8,24]  with  a  Hamiltonian  connecting  quantum  rotors 
on  a  lattice  with  /3S  acting  as  the  coupling  between  the  spatial 


sites.  In  the  ®A  \nx)  basis,  it  reads 

cosWx  -  0y),  (3) 

x  x  ( xy ) 

with  0  =  1  / (ySr «),  /x.  =  ii/a,  and  J  =  ps/a,  the  sum  extend¬ 
ing  over  sites  x  and  nearest  neighbors  (xy)  of  the  space  lattice, 
and  a  is  a  lattice  spacing. 

The  commutation  relations  [ L,e±w \  —  ±e±,e  show  that 
e±,e  act  like  creation  and  annihilation  operators.  However, 
there  is  no  eigenstate  of  L  annihilated  by  e~'e .  At  large  /z,  there 
is  an  effective  truncation  [25,26]  that  makes  the  eigenstates 
with  negative  eigenvalues  irrelevant.  For  small  values  of  /x, 
we  will  consider  the  quantum  link  inspired  truncation  where 
the  original  operator  algebra  is  approximated  by  a  spin-.s 
representation  with  \n\  ^  s. 

Remembering  the  role  played  by  the  differential  operator 
L  =  —id/d9  in  the  construction  of  the  spherical  harmonics, 
we  replace  L  by  iJ ,  the  third  component  of  the  angular 
momentum  in  the  SU(2)  Lie  algebra.  Pursuing  the  analogy, 
we  replace  e±,e  by  an  operator  proportional  to  the  raising 
and  lowering  operators  L*  in  the  spin-,y  representation.  In  the 
case  of  spin-1,  a  comparison  of  the  matrix  elements  shows 
that  the  correspondence  between  the  two  representations 
can  be  accomplished  by  properly  choosing  the  constant  of 
proportionality. 

HI.  NUMERICAL  CALCULATION  OF  THE 

PHASE  DIAGRAM 

We  now  discuss  the  phase  diagram,  the  finite  spin  projec¬ 
tion,  and  the  time  continuum  limit  by  using  the  TRG  method. 
Following  the  procedure  described  in  Refs.  [8-10],  we  can 
write 

Z  =  Trn  T(x'?  , ,  (4) 

1  I  nxnxntnt  ’  v  / 

(X.t) 

with  the  local  tensor  expressed  in  terms  of  the  modified  Bessel 
functions: 

Tn*n%tnt,  =  exp[^(n,  +  <)] 

X  yj ' Itix  (fis)^nx+nt ,nx>+nti  •  (5) 

The  indices  nx,  n'x.  n,.  and  n't  label  the  four  links  coming  out 
of  (x,t)  in  the  x  and  t  direction,  and  the  trace  Tr  refers  to  the 
sum  over  all  these  link  indices.  A  transfer  matrix  T  can  be 
constructed  by  taking  the  spatial  traces  in  a  time  slice: 

TT  ill  _  \  '  'T’(lU)  t42,0 

^  1  >ixNsnx\n\ni  nxlnx2n2n2*--- 


The  indices  (n\ .ni,  . . .  ,»jvs)  represent  the  past  and 
(n\  ,n'2,  . . .  ,n'N  )  the  future. 

In  view  of  the  rapid  decay  of  the  I„(fi)  when  \n\  increases  at 
fixed  /3,  good  approximations  can  be  obtained  by  replacing  the 
infinite  sums  by  sums  restricted  to  —  nmax  to  ;imax.  We  denote 
the  number  of  states  Dst  =  2nmax  +  1.  With  this  truncation, 
the  transfer  matrix  is  a  D[//  x  D\//  matrix.  It  is  possible  to 
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TABLE  I.  (A)  for  the  worm  algorithm  and  the  HOTRG  for  f)s  = 

Pr  =P- 


fi 

M 

(A)  (worm) 

(A)  (HOTRG) 

1.12 

0.01 

0.00726(1) 

0.00728(8) 

0.46 

1.8 

0.98929(1) 

0.9892(3) 

0.28 

2.85 

1.98980(2) 

1.989(2) 

0.2 

3.53 

2.96646(3) 

2.967(1) 

0.12 

4.3 

3.96206(4) 

3.965(1) 

coarse  grain  the  transfer  matrix  efficiently  by  using  a  higher- 
order  singular  value  decomposition  (HOTRG)  described  in 
Ref.  [7].  This  procedure  then  reduces  the  two-site  transfer 
matrix  to  a  Dst  x  Dst  matrix  and  thus  accomplishes  the 
blocking  from  two  sites  to  a  single  site.  Note  that  in  the  spin- 1 
projection  we  keep  Ds,  much  larger  than  3  as  we  keep  blocking. 
In  other  words,  the  spin  projection  represents  a  microscopic 
modification  of  the  model,  while  we  need  to  keep  Dst  as  large 
as  possible  in  order  to  maintain  good  macroscopic  accuracy. 
The  same  numerical  method  is  used  in  all  cases,  the  only 
difference  being  the  initial  tensor. 

An  important  advantage  of  the  TRG  method  is  that  it 
allows  us  to  reach  exponentially  large  volumes.  However,  it  is 
important  to  check  the  results  at  small  volume  where  sampling 
methods  are  feasible  and  accurate.  We  have  used  the  TRG  and 
the  worm  algorithm  [11,12]  to  calculate  the  particle  number 
density  [12] 

(N)  =  1/(77,  x  Nz)d  In  Z/9/x.  (7) 

The  partition  function  Z  can  be  calculated  by  taking  the  trace 
of  or  by  using  the  methods  described  in  Refs.  [7-10]. 
The  numerical  values  on  a  16  x  16  lattice,  for  values  of  fis  = 
fiT  and  jjL  slightly  below  the  tips  of  the  regions  with  (N)  = 
0, 1, ...  ,4  of  the  phase  diagram  described  below,  are  shown  in 
Table  I. 

Small  discrepancies  between  the  two  methods  appear 
typically  in  the  fourth  significant  digit.  The  errors  for  the 
worm  algorithm  are  purely  statistical,  and  to  the  best  of  our 
knowledge,  there  are  no  systematic  errors  associated  with  it. 
On  the  other  hand,  for  the  TRG  method,  the  limit  Dst  oo 
shows  very  small  variations,  which  will  be  documented  and 
analyzed  in  a  separate  publication  [27]  but  do  not  affect  the 
results  presented  here. 

By  increasing  \x  at  fixed  /3,  we  go  through  successive 
Mott  insulating  (MI)  phases  characterized  by  a  fixed  integer 
value  of  (A)  increasing  with  /x  and  alternating  with  superfluid 
(SF)  phases  where  (A)  interpolates  continuously  between  the 
successive  integers.  The  phase  boundaries  are  clearly  visible 
from  the  steps  in  ( A )  as  a  function  of  /x,  as  shown  in  Fig.  1. 
The  phase  boundaries  can  also  be  obtained  by  looking  at  the 
two  largest  eigenvalues  of  the  transfer  matrix.  In  a  given  MI 
phase,  one  would  expect  that  the  largest  value  of  the  transfer 
matrix  is  unique  and  corresponds  to  an  eigenstate  with  fixed 
integer  particle  density.  On  the  other  hand,  in  the  SF  phase,  the 
two  largest  eigenvalues  of  the  transfer  matrix  are  expected  to 
be  degenerate  and  the  corresponding  eigenstates  are  expected 
to  have  particle  density  corresponding  to  the  two  neighboring 
MI  regions.  Figure  1  shows  that  these  expectations  are  verified 
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B 

FIG.  1 .  (Color  online)  The  ratio  (At  /ki)  of  the  first  two  eigenval¬ 
ues  of  the  transfer  matrix  and  the  particle  number  density  (A)  for  the 
fiT  =  //  =  0.06  from  the  HOTRG  calculation  with  Dst  =  15.  The 
particle  number  density  (A)3  and  the  second  normalized  eigenvalues 
(k2/k  1)3,  where  a  lower  index  3  denotes  the  spin-1  projection 
(three-state),  are  also  shown. 

quite  precisely.  The  system  reaches  the  superfluid  (SF)  phase 
when  A.2/A  1  =  1,  and  when  /x  is  increased  further,  this  ratio 
remains  1  while  there  is  an  increase  in  the  particle  number 
density  between  two  adjacent  integers,  which  stand  for  two 
different  MI  phases. 

The  alternation  between  the  MI  and  SF  phases  in  the 
/3-/x  plane  is  shown  in  Fig.  2.  The  pointy  shape  of  the  MI 
phase  region  is  also  observed  in  other  (1  +  l)-dimensional 
Bose-Hubbard  models  [28,29].  The  spin-1  projection  is  also 
shown  in  these  figures.  When  /x  is  not  too  large,  only  small 


FIG.  2.  (Color  online)  The  phase  diagram  in  the  f) -/x  plane  for 
the  isotropic  (/?,  =  f)T  =  /3)  case. 
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differences  with  the  original,  unprojected  model  are  observed. 
However,  when  //  becomes  large  enough  to  have  (N)  >  1,  the 
truncation  prevents  such  a  large  occupation  and  (N)  saturates 
to  1  as  expected,  and  there  is  no  (N)  =  2  MI  phase.  The  phase 
boundary  in  Fig.  2  between  the  MI  (N)  =  0  phase  from  the  SF 
phase  coincides  approximately  with  the  line  for  the  model  with 
an  infinite  number  of  states.  Similarly,  the  spin-2  projection 
(not  shown  on  the  graph)  reproduces  well  the  (N)  =  0  and  1 
boundaries  while  discrepancies  appear  for  (N)  =  2. 

Figure  2  shows  that  when  p  is  small,  the  boundary  between 
the  MI  and  SF  phase  appears  to  be  at  large  values  of  /x.  It  is 
useful  to  recall  that  so  far  we  have  only  considered  the  phase 
diagram  in  the  isotropic  case  ft  =  ps  =  /3Z.  When  p>  — >  0, 
the  interaction  along  the  space  links  is  small,  but  if  /x  is 
sufficiently  large,  the  interactions  along  the  time  links  are  not 
small.  In  the  limit  where  the  interactions  among  the  space 
links  are  negligible,  the  problem  reduces  to  a  collection  of 
independent  one-site  problems  (simple  quantum  mechanics) 
as  in  mean-field  theory  [14].  In  this  limit,  Eq.  (5)  shows  that 
the  transfer  matrix  becomes  diagonal  because  7„(0)  =  0  except 
for  n  =  0  [7o(0)  =  1],  and  by  the  conservation  law  the  same 
index  nx  characterizes  the  interaction  along  the  time  direction. 
In  other  words,  there  is  no  quantum  number  flowing  in  the 
space  direction,  and  the  flow  in  the  time  direction  at  each  site 
is  constant.  In  this  limit,  the  eigenvalues  of  the  transfer  matrix 
are  just 

hnl,n2,...,nNs)  =  l\lnx(P)en*11.  (8) 

X 

The  largest  eigenvalue  is  then 

kmax  =  {ma xn[In(P)enn}Ns  ■  (9) 

Finding  the  value  of  n  corresponding  to  the  maximum 
eigenvalue  gives  the  particle  density  (N)  in  the  MI  phase. 

The  maximization  of  An  =  l„(p)enl‘  can  be  achieved  by 
considering  the  ratios  rn  =  An+\/An.  Note  that  we  assume 
/x  >  0,  and  given  that  I„(P)  =  I-n(P),  we  only  need  to 
consider  n  p  0.  Whenr„_i  >  1  and  r„  <  1 ,  A, ,  is  a  maximum. 
It  can  be  shown  in  the  limit  of  small  and  large  p  that  rn 
decreases  when  n  increases.  If  this  remains  true  for  arbitrary 
P  and  if  rn  ^  1,  then  the  problem  has  a  unique  solution. 
The  interesting  case  is  rn  =  1,  which  implies  An  =  An+\  and 
should  be  at  the  boundary  between  two  MI  phases  with  particle 
density  n  and  n  +  1.  In  the  small-/?  limit,  I„(P)  —  pn/(2"n\) 
and  the  condition  r„  =  1  implies  that 

Pe^t  2  =  n+  1  (10) 

in  that  approximation.  The  sudden  transition  in  particle  density 
occurs  at  integer  values  of  /3<?M/2.  This  prediction  is  confirmed 
by  plotting  the  phase  diagram  in  the  /?-/3eM/2  plane  as  shown 
in  Fig.  3.  We  see  that  by  changing  the  vertical  coordinate  to 
/x  — »■  PeJL  / 2,  the  shape  of  the  phase  diagram  of  the  isotropic 
system  looks  like  the  cuspy  shapes  found  for  the  Bose-Hubbard 
model  in  one  spatial  dimension  [28,29],  Keeping  in  mind  that 
we  are  working  in  the  limit  of  small  P,  Eq.  (10)  implies  that 
the  phase  boundaries  of  the  SF  phase  between  the  n  and  n  +  1 
MI  phases  diverge  like  ln[2(«  +  1 )//?]  when  p  -+  0,  which  is 
consistent  with  Fig.  2. 

We  now  depart  from  the  isotropic  pT  =  ps  situation 
and  consider  the  case  pT  'A>  ps  corresponding  to  the  time 


FIG.  3.  (Color  online)  The  phase  diagram  in  the  /6-J6efl/2  plane 
for  the  isotropic  case.  ( N )  =  0  (/x  =  0)  line  in  the  SF  phase  is  the 
Kosterlitz-Thouless  phase  in  the  1  +  ID  0(2)  model  at  fi  =  0.  The 
lines  labeled  by  “3s”  stand  for  the  phase  boundaries  of  the  spin-1 
(three-state)  system. 

continuum  limit.  If  we  neglect  ps,  we  obtain  the  one-site 
approximation  described  above.  The  particle  density  can  be 
obtained  from  the  ratio  analysis  in  the  larg e-/?T  limit.  Using 
In+i(PT)/ In(PT)  —  1  —  (n  +  1/2 )/pz  in  this  limit,  we  find 
that  the  degeneracy  occurs  for  integer  values  of  /xyST  —  1/2. 
Defining  the  effective  chemical  potential  /xf  =  iipT  —  1/2  and 
effective  coupling  pe  =  psPi,  we  find  that  the  same  MI-SF 
pattern  appears  in  the  pe-/ie  plane  (Fig.  4). 

Having  computed  the  phase  diagram  in  the  /?T  =  ps  and 
Pz  2S>  Ps  cases,  we  learned  that  they  have  very  similar  shapes 
in  suitable  systems  of  coordinates.  From  this  we  expect  that 


PtP: 


FIG.  4.  (Color  online)  The  phase  diagram  for  the  1  +  ID  0(2) 
model  at  fir  =  10  in  the  pe-^ie  plane  is  shown.  The  lines  labeled  by 
“3s”  stand  for  the  phase  boundaries  of  the  spin- 1  (three-state)  system. 
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they  can  be  smoothly  deformed  into  each  other  and  that  nothing 
special  happens  in  the  intermediate  situations. 

IV.  OPTICAL  LATTICE  IMPLEMENTATION 

To  incorporate  the  positive  and  negative  eigenvalues  of  L, 
we  will  consider  a  two-species  Bose-Hubbard  Hamiltonian  on 
a  lattice: 

H  =  -  YM*  +  tbb*by  +  Hc)  “  +  Aa)n“ 

( xy )  x,oc 

+J2  -  0  +  w^2n>bx  +  J2 

x,ce  x  (xy),a 

(ii) 

with  a  =  a,b  indicating  the  two  different  species,  nax  —  a\ax 
and  nbx  =  b\bx  the  number  operators,  and  \nx,nb)  the  corre¬ 
sponding  on-site  basis.  This  class  of  models  has  been  studied 
extensively  [30-33].  It  is  possible  to  adjust  the  chemical 
potentials  in  order  to  set  ( nx )  =  (nx  +  nb)  =  2.  In  the  limit 
where  Ua  =  U\,  =  W  are  very  large  and  positive,  the  on-site 
Hilbert  space  can  then  be  restricted  to  the  states  satisfying 
nx  =  2  at  each  site.  All  the  other  states  (with  nx  ^  2)  belong  to 
high-energy  sectors  that  are  separated  from  this  one  by  energies 
of  order  U.  The  three  states  |2,0),  1 1, 1),  and  |0,2)  correspond 
to  the  three  states  of  the  spin-1  projection  considered  above. 

It  is  useful  to  visualize  the  minima  of  the  on-site  Hamilto¬ 
nian  obtained  in  the  limit  t  —>  0.  It  can  be  written  as  a  quadratic 
form  and  a  linear  term  in  na  and  «/,.  If  UaUb  >  W2,  there  is 
a  unique  minimum,  1 1, 1),  which  corresponds  to  a  miscible 
phase  where  the  two  species  need  to  be  present  at  the  same 
place.  Since  in  the  spin-1  approximation,  1 1 , 1)  corresponds 
to  a  rotor  with  angular  momentum  zero,  this  is  the  correct 
situation  for  the  0(2)  model  we  try  to  simulate.  On  the  other 
hand,  if  UaUb  <  VI'2.  the  extremum  is  a  saddle  point.  As  we 
will  discuss  later,  the  unstable  direction  coming  out  of  the 
extremum  is  limited  by  the  positivity  of  the  occupation  number. 
There  are  two  vacua  |2,0)  and  |0,2),  which  corresponds  to 
immiscible  phases.  The  limiting  case  UaUb  —  W  corresponds 
to  our  Ua  =  Ub  =  W  =  Uq  lowest-order  approximation.  If  in 
addition  we  have  ji  =  (3/2)t/o  and  Aa  =  0,  we  have  a  flat 
direction  along  the  line  nx  =  2  where  we  have  three  states 
of  energy  —2Uq,  while  the  degenerate  lines  with  nx  =  I  or 
3  have  energy  —  (3/2)£/o,  which  is  considered  much  larger 
in  the  strong-coupling  approximation.  Small  changes  in  the 
parameters  will  break  the  degeneracy  of  the  ground  state  but 
preserve  a  significant  difference  between  these  states  and  the 
excited  states.  Decreasing  W  lowers  the  energy  of  the  1 1 , 1) 
state  linearly  in  the  difference  with  Uo ■  Similarly,  increasing 
W  raises  the  energy  of  the  1 1 , 1 ) ,  the  flat  direction  curves 
down  at  both  ends,  but  the  positivity  of  the  occupation  number 
prevents  the  energy  from  being  unbounded  from  below.  When 
species-dependent  chemical  potentials  are  turned  on,  the  flat 
direction  becomes  slanted  linearly  in  the  variation  of  the 
chemical  potential  Aa.  The  overall  shape  of  the  trap  will 
typically  create  small  variations  in  a  space-dependent  manner. 
In  summary,  as  long  as  the  variations  of  the  parameters 
are  small  compared  to  Uo,  the  features  departing  from  the 
degenerate  case  can  be  treated  as  perturbations. 
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FIG.  5.  (Color  online)  Two  species  (green  and  red)  of  bosons  on 
species-dependent  optical  lattices  (with  the  same  color).  The  nearest- 
neighbor  interaction  is  coming  from  the  overlap  of  Wannier  Gaussian 
wave  functions.  We  assume  the  difference  between  inhaspecies 
interactions  are  small  U  8. 

Going  back  to  the  general  Hamiltonian  [Eq.  (1 1)],  we  write 
Ua(b)  =  U  ±  8  and  assume  U  S,(U  —  W),V  ,ta,  Aa  and  do 
degenerate  perturbation  theory.  Virtual  processes  exchanging 
particles  between  neighboring  sites  are  allowed  at  second  order 
with  contributions  proportional  to  —tata'/U.  The  hopping 
amplitude  is  tunable,  and  when  chosen  to  be  ta  =  *jVaU /2, 
the  final  result  is  that  the  effective  Hamiltonian  up  to  second 
order  in  degenerate  perturbation  theory  corresponds  to  the 
spin-1  projection  of  the  rotor  Hamiltonian  of  Eq.  (3)  with 
J  =  4YaVb,  0  =  2(U  —  W),  and  /l  =  (Aa  -  Va)  -  (A*  - 
Vb).  Similarly,  by  increasing  the  chemical  potentials,  it  is 
possible  to  restrict  the  Hilbert  space  to  nx  +  nb  —  2s,  which 
corresponds  to  a  spin-.v  projection  in  the  0(2)  model. 

This  two-species  Bose-Hubbard  model  can  be  realized  in 
a  87 Rb  and  41 K  Bose-Bose  mixture  where  an  interspecies 
Feshbach  resonance  is  accessible  [34,35],  Due  to  the  physical 
nature  of  the  different  atoms,  the  hopping  amplitudes  (ta ,  f/,)  are 
different  to  begin  with,  as  well  as  the  intraspecies  interactions. 
In  addition,  species-dependent  optical  lattices  [36-40]  are 
widely  used  in  boson  systems,  which  allows  the  hopping 
amplitudes  of  each  individual  species  to  be  further  tuned  to  the 
desired  value.  As  mentioned  above,  the  interspecies  interaction 
(W)  can  be  controlled  by  an  external  magnetic  field  [35]. 
Finally,  the  extended  repulsion,  Va,  is  present  and  small  when 
we  consider  Wannier  Gaussian  wave  functions  centered  on 
nearby  lattice  sites  according  to  previous  study  [41].  This 
is  schematically  illustrated  in  Fig.  5.  This  may  be  the  most 
difficult  parameter  to  achieve,  but  other  proposals  may  be 
explored,  such  as  by  using  dipolar  bosons  [42],  or  by  pumping 
bosons  to  higher  Bloch  bands  [43]  in  order  to  engineer  the 
nearest-neighbor  interaction.  It  is  also  important  to  have  U 
significantly  larger  than  the  temperature.  For  the  mixture 
considered  here,  the  temperature  and  recoil  energies  are  of 
the  order  of  100  nK,  and  values  of  U  10-20  times  larger  can 
typically  be  reached  [35,44,45]. 

V.  CONCLUSIONS 

In  summary,  we  have  used  recently  developed  numerical 
methods  to  connect  the  0(2)  model  in  1  +  1  dimensions  to 
an  optical  lattice  setup.  A  first  test  of  the  correspondence 
would  be  to  check  that  the  optical  lattice  system  reproduces 
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the  phase  diagram  of  Fig.  4,  which  corresponds  to  the  time 
continuum  limit  (iT  fts  of  the  classical  model  and  where 
the  microscopic  parameters  can  be  approximately  connected 
to  those  of  the  two-species  Hubbard  model. 

The  TRG  method  presented  here  allows  reliable  calcula¬ 
tions  of  the  eigenvalues  X,  of  the  transfer  matrix.  In  the  time 
continuum  limit,  we  have 

Xi/Xl  oc  e-a(E‘~E°\  (12) 

with  E,  the  corresponding  energies  and  a  oc  1  / fiT  the  lattice 
spacing.  Recently  developed  experimental  techniques,  e.g., 
momentum-resolved  Bragg  spectroscopy  [46],  could  in  prin¬ 
ciple  allow  detailed  comparisons. 

We  have  shown  that  for  low  enough  //,  the  effect  of  the 
truncation  to  spin-1  or  -2  of  the  original  0(2)  model  had  small 
effects  on  the  phase  boundaries.  In  the  TRG  formulation,  this 
truncation  only  affects  the  initial  values  of  the  tensor,  which 
can  be  compared  with  the  initial  tensor  of  other  spin  models 
with  a  finite  number  of  states  in  configuration  space  (clock 
and  Potts  models).  Understanding  how  the  symmetries  of  this 
initial  tensor  affect  the  universality  class  is  under  study. 

The  0(2)  model  has  an  exact  conservation  law  that  is  made 
clear  by  the  Kronecker  <5  in  Eq.  (5).  The  states  kept  in  the  TRG 
calculation  have  a  well-defined  quantum  number  associated 
with  this  conservation  law,  and  it  can  monitored  and  put  into 
histograms  [27],  This  provides  detailed  information  about  the 


average  occupation  and  its  fluctuations.  It  could  give  better 
insight  into  the  validity  of  the  Gutzwiller  ansatz  or  mean-field 
calculations  such  as  the  ones  discussed  in  Ref.  [14],  or  the 
validity  of  the  finite  spin  projection  discussed  here. 

In  LGT  calculations,  important  information  regarding  the 
spectrum  and  matrix  elements  can  be  extracted  from  the  two- 
and  three-point  functions  obtained  by  introducing  localized 
sources  in  the  Lagrangian  formulations.  Techniques  to  gather 
related  information  from  an  optical  lattice  system  have  yet  to 
be  developed.  Generalizing  the  work  done  here  for  the  0(3) 
model,  which  has  physics  more  similar  to  lattice  QCD,  seems 
possible  and  interesting. 
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